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REAL-TIME SIMULTANEOUS ANALYSIS AND SHORT ti=dm 
PREDICTION OF HEART AND PULSE^ATE VARIAB^LnT 

FIELD OF THE INVENTION 

The present invention relates to the provision of non-stationary and 
non-linear heart rate variability data.and more particularlyto immediate processing 
analys.s and display of heart rate variability data and its surrogates using non- 
invasive analysis techniques. 

GOVERNMENT SUPPORT 

The present invention was made with U.S. Government support from 
the National Institutes of Health, National Heart, Lung, and Blood Institute, under 
Grant No. HL 67735. and the National Institute of Neurological Diseases and 
Stroke, under Grant No. NS 37981 . The U.S. Government has certain rights in this 
invention. 

BACKGROUND 

Measurements of heart rate and its variability are well known in the art 
for their usefulness in assessing the conditions of the circulatory system in both 
health and in disease. They are useful for monitoring many chronic diseases, such 
as diabetes and heart failure, as well as for monitoring cardiac status during exer- 
cise. Particularly useful is Heart Rate Variability (HRV) analysis, which is a non- 
invasive, clinical tool for assessing the autonomic regulation of cardiac activity. It is 
also used as a marker of psychophysiological function for biofeedback therapy. 
Reduced HRV has been associated with such problems as higher long-term risk of 
post-infarction mortality. 
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HRV analysis is commonly performed by measuring the beat-to-beat 
interval between successive heartbeats as collected on an electrocardiogram 
(ECG). A particularly useful parameter is the period between succeeding "R" 
waves (the R-R interval), where "R" is the conventional designation given the 
5 waveform peak of a normal heartbeat 1 0 as illustrated in Fig. 1 . Most analyses of 
short-term electrocardiograms using conventional frequency domain HRV tech- 
niques (e.g., power spectral density) assume "stationarity" of the underlying R-R 
interval time series. However, most physiological signals, including heart rate (HR) 
and pulse rate (PR), are non-stationary by nature. This non-stationarity is a result 
1 0 of complex dynamic interactions among multiple bioregulatory control mechanisms 
responsible for maintaining homeostasis in the presence of constantly varying 
physiological and environmental inputs. Additionally, conventional spectral analysis 
methods are limited by their inability to assess transient changes in HR and PR 
associated with temporary physical or mental stresses or cardiac pathologies. 
1 5 Joint time-frequency (t-f) signal processing techniques may be advan- 

tageously used over conventional tools for HRV analysis given their ability to ana- 
lyze time-varying spectral properties of non-stationary signals such as HRV. Such 
t-f techniques are ideally suited for time-localized spectral characteristics of tran- 
sient cardiac events that occur as a result of temporal changes in the sympatho- 
20 vagal balance. However, the common use of the Gabor spectrogram, where the 
short-time Fourier transform is calculated for a window chosen to be appropriate for 
the data to be collected, may make it difficult to achieve an appropriate compro- 
mise between frequency resolution and time resolution. 

Techniques such as chaotic analysis have the ability to assess non- 
25 linear, spatio-temporal behavior of such deterministic systems as cardiac activity. 
Additionally, chaotic analysis has the potential for predictive value in the screening 
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of patients susceptible to lethal arrhythmias. A numerical "chaotic index" derived 
from the non-linear analysis of an electrocardiogram or a pulse-rate signal can be 
used to quantify the degree of non-linear deterministic behavior of cardiac activity. 
Techniques developed out of chaos theory, such as embedding methods and 
estimation of Lyapunov exponents, help to unravel the original signal underlying an 
observed single-variable time series and determine how far into the future it can be 
predicted. Chaotic systems comprise a class of signals that lies between predict- 
able periodic or quasi-periodic signals and totally irregular stochastic signals which 
are completely unpredictable. The Lyapunov exponent measures the sensitivity of 
the system to initial conditions and thus provides a measure to help predict the 
short-term behavior of the system. 

Such HRV analysis has heretofore typically been performed by hook- 
ing a patient up to equipment for monitoring heart activity and storing the data from 
the monitored heart activity. After monitoring heart activity has been accomplished 
for a sufficient period of time so that a selected amount of data has been accumu- 
lated (e.g., from several minutes to several hours), the data are transferred to a 
computer in which they are analyzed to provide the physician with such information 
as the Chaotic Index (the largest Lyapunov exponent [measure of degree of chaos] 
calculated using the data represented by the heart rate [or pulse rate] sequence), 
BPM (beats per minute [Heart Rate or Pulse Rate]), SDNN (standard deviation of 
R-R intervals [or inter-pulse intervals] derived from the electrocardiogram [or pulse 
plethysmograph] data after putative abnormal R-R intervals [or inter-pulse intervals] 
are removed), and RMSSD (root-mean-square of the difference between succes- 
sive R-R intervals [or inter-pulse intervals] derived from the electrocardiogram [or 
pulse plethysmograph] data). The generated information is then reviewed by a 
physician, typically long after the heart activity which was used to generate the 
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information has taken place, and the physician uses the generated information at 
that later time to determine a status or treatment procedure for the patient. Task 
Force of the European Society of Cardiology and the North American Society of 
Pacing and Electrophysiology. Heart Rate Variability: Standards of Measurement, 
5 Physiological Interpretation, and Clinical Use. Circulation, 93(5), pp1 043-1 065, 
1996; Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh. Mark RG, 
Mietus JE, Moody GB, Peng CK, Stanley HE. PhysioBank, PhysioToolkit, and 
Physionet: Components of a New Research Resource for Complex Physiologic 
Signals. Circulation 101(23): e215-e220 and U.S. Patent Nos. 5,265,617, 

10 5,437,285, 5,682,901 , 5,842,997, 5,957,855, 6,1 1 5,629, 6,41 6.471 , 6,480,733, and 
6,485,416 variously teach HR monitoring and analysis, and their full disclosures are 
hereby incorporated by reference. 

While such analysis has had value in the treatment of a patient, the 
delay in the analyzed data provided to the physician has clear disadvantages. For 

1 5 example, the physician's receipt of analyzed data may be so delayed as to cause 
the initiation of an action, such as a treatment, to be disadvantageous^ delayed, or 
in the worst case the information may be generated too late to be of help in treating 
the patient. As another example, the review of such information by a physician 
hours after the data were collected may make it difficult to correlate the data with 

20 other conditions of the patient for which data were not being simultaneously re- 
corded. A physician will observe many things when with a patient, but hours later 
may be unable to temporally correlate many of those observations with the corre- 
sponding HRV data. Still further, it should be appreciated that prior art HRV infer, . . . 
mation which has been generated based on a selected set of data has in many 

25 ways been able to present only a static picture of a dynamic situation. 
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The present invention is directed toward overcoming one or more of 
the problems set forth above. 



SUMMARY OF THE INVENTION 

Since HRV is an important measure of the condition of the heart, 
5 which could change rapidly, it is particularly important to be able to perform real- 
time HRV analysis while the heart rate is being measured. No current system is 
able to perform such detailed HRV analysis and display the results in real time. 
The present invention relates to a system and method that can simultaneously 
acquire electrocardiogram or pulse rate data, dynamically perform time-frequency 

10 (t-f) and chaotic analysis in real-time, visually display the results in a convenient 
graphical format and store the results in a computer file format. The system and 
method can provide a real-time automated system that combines the non-station- 
ary analysis capability for evaluating cardiac signal histories with a predictive capa- 
bility of non-linear analysis to better monitor and categorize autonomic regulation 

15 of cardiac function. This system allows for continuous, real-time monitoring of 
cardiac function and enables short-term prediction of cardiac activity under auto- 
nomic control. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Figure 1 illustrates a waveform of a conventional heart beat; 
20 Figure 2 is an example computer and data acquisition device which 

may be used in accordance with the present invention; 

Figure 3 is a monitor showing an example screenshbt from the sys- 
tem performing HRV analysis on pre-acquired representative data for a normal 
electrocardiogram; 
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Figure 4 is a flowchart illustrating the overall process of data acquisi- 
tion, analysis and display of the real-time HRV or pulse rate variability (PRV) analy- 
sis system according to the present invention, where Figs. 5-9 are detailed flow- 
charts of portions of the Fig. 4 process: 
5 Figure 5 is a flowchart illustrating the initiation of analysis and data acquisi- 

tion; 

Figure 6 is a flowchart illustrating details of the event detection step; 
Figure 7 is a flowchart illustrating details of heart rate resampling and RR 
sequence generation; 

10 Figure 8 is a flowchart illustrating details of determining time-frequency 

distribution; and 

Figure 9 is a flowchart illustrating details of non-linear data analysis; 

Figures 10A-10C are sequential pages of an algorithm in Visual C++ 
programming language which may be used to detect the peaks in the R wave and 
15 to compute the R-R interval; . 

Figure 1 1 is an example screenshot from the system performing HRV 
analysis on pre-acquired representative data for a Chronic Heart Failure (CHF) 
electrocardiogram; and 

Figure 12 is an example screenshot from the system performing HRV 
20 analysis on pre-acquired data representative of an epileptic seizure episode elec- 
trocardiogram. 

The accompanying drawings, which are incorporated in and form a 
part of this specification, illustrate embodiments of the invention and, together with 
descriptions, serve to explain the principles of the invention. They are not intended 
25 to limit the scope of the invention to the embodiments described- It is appreciated 
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that various changes and modifications can be made without departing from the 
spirit and scope of the invention as defined in the appended claims. 

DETAILED DESCRIPTION OF THE INVENTION 

Fig. 2 illustrates a device which may be used in accordance with the 
5 present invention. As illustrated, the device comprises a suitable processing unit, 
such as a personal computer 12 with a suitable CPU, a user input device 14 (such 
as the illustrated keyboard, and/or other suitable devices such as a mouse, touch- 
screen, or keypad-controlled graphic user interface)* and suitable display device 
such as a CRT monitor 16, and a suitable data acquisition device 18 which may be 

1 0 attached to a patient to obtain ECG data of a patient's heart. 

In order to facilitate processing of the various elements of HR V analy- 
sis which may be performed in real-time in accordance with the present invention, 
a suitable processor with architecture for performing specific functions may be 
advantageously used with the computer 12. For example, Intel Corporation's IPP 

15 ("Integrated Performance Primitives") for its Pentium® processors and Itanium® 
architecture permit a variety of operations which are performed in connection with 
the present invention to be quickly performed, and may therefore be advanta- 
geously used in a computer 12 used with the present invention. For example, 
where Visual C++ based software code is used to perform the operations providing 

20 the desired HRV analysis, the following operations may be performed using IPP 
function calls: memory allocation and deallocation, array initialization, freeing of 
memory, calculating means, absolute values and exponentials for array elements 
in multidimensional arrays, Fast Fourier Transforms and Inverse Fast Fourier 
Transforms. Use of such IPP function calls permit such functions to be performed 

25 with significantly fewer lines of software code than would be required to perform 
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such functions with normal processing, and therefore can significantly speed up the 
processing functions to allow analysis to stay in real-time as the processes are 
continuously updated as discussed in further detail below. 

Suitable ECG data acquisition devices 18 acquire heart beat data 
5 such as is known in the art, and are available from, for example, QRS Diagnostic, 
LLC of Plymouth, Minnesota, U.S.A., which have devices which may be connected 
directly to a serial port on a PC processing unit without requiring special hardware 
to communicate through the serial port. However, it should be understood that 
many different data acquisition devices may also be advantageously used within 

10 the scope of the Invention including, for example, devices which may be directly 
connected to different computer ports, such as PCMCIA ports conventionally found 
in laptop computers. In addition, a battery-powered device, such as a cell phone, 
PDA, or tablet PC may be used in combination with an electrocardiogram or finger 
pulse sensor, for acquisition, transmission and remote monitoring of heart rate 

1 5 variability or pulse rate variability parameters. Further, it should be recognized that 
multiple devices may be used with a single computer (whether other computer 
components or other ECG data acquisition devices), with the connection to the 
ECG data acquisition device of interest (connected to the patient of interest) being 
selectable by the user. As another example, it should be appreciated that the 

20 present invention may utilize a battery-powered, Class H biofeedback prescription 
device, such as a PDA or a battery powered computer, in combination with a finger 
pulse sensor or ECG acquisition system. 

Further, as detailed herein, the data acquisition device 18 may be 
used in connection with the present invention to provide data for real-time analysis 

25 simultaneously with its collection, or may acquire data which is suitably stored (e.g., 
on the hard drive of a personal computer 1 2) in a form which preserves moment-to- 
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moment correlative relationships so that they can be retrieved electronically for 
later playback or as hard copy for later review and/or documentation. 

Fig. 3 is an example of a video display on a monitor 16 of the HRV 
analysis of a normal ECG in accordance with the present invention, where the data 
5 being analyzed have been pre-acquired: 

Graphic display element 20 illustrates in wave form the unprocessed elec- 
trocardiogram or pulse data of the pre-acquired data file; 
Graphic display element 22 is the heart rate or pulse rate derived from the 
electrocardiogram or pulse data; 
1 0 Graphic display element 24 is an intensity-mapped time-frequency distribu- 

tion color contour plot, with its time axis (the horizontal axis) shared 
with the time axis of display element 22; 
Graphic display element 26 is the electrocardiogram or pulse attractor 
("ECG attractor") derived from the electrocardiogram or pulse data; 
1 5 Graphic display element 28 is the play button for pre-acquired data process- 

ing (if operated while the data are being acquired from a patient, the 
element 28 functions as a start button); 
Graphic display element 30 is the stop button both for pre-acquired data 
processing and for the mode of operation in which the data are pro- 
20 cessed real-time as it is acquired from a patient; and 

Graphic display element 32 is the dialog box for reporting system status. 
In the example shown, the heart rate is 1 17 beats per minute, the SDNN is 88 ms, 
the RMSSD is 22 ms, the SI/PI ratio (discussed further below) is 1.3, and the 
Chaotic Index is 0.41. Moreover, as detailed hereafter, this information is dynamic 
25 in that it is analyzed and modified over real time. That is, even when used with a 
file of pre-acquired data, the analysis and display of the data and analysis occurs 
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dynamically over „me corresponding «o the passage of flme which ocourred when 
tee date were acquired. ,« should a,so be appreciated that when used as the data 
« -being acquired, as discussed herein. tee analysis and dispiay o, the Z « 
ana,ys,s. occurs dynamicaiiy as tee dynamic event (,.e, patten, hear, heat, occure 
F,g. 4 ,s a flow chart illustrating tee dynamic, real time HRV analysis 
wh,chmay b epertormedinacco ra ancew«htee pre5 en.i„venflon. Further^" 
of te« opemhon ...ustreted in overview in Fig. 4 are se, forth hereinafter, including 
in Figs. 5-9 and in the associated written specification. 

detailed, rth ^T^" " ^ 40 *" ^ '*» P8R " "* « 
detailed further befow. Such date parametere may be used to cent™, tee analyse 

mode and output generation, including whether analysis is to be performed using 

ane«ema.dataseurce(a t 4 2 )orafl,eofp re .ac q uirede,ec,roca ra iogramerpul 8 e 
waveform date (a. 44). Such parametere can further include, for example, sam- 
pling frequency for real time date acquisition, or selection of the file which has the 
pre-acqu,red date of interest. The date to be used (such as unprocessed ECG 
date or pu.se plethysmograph data) and the parameters to use in connection with 
its ana.ysis are teen received at box 46. This process is set forth in greater detail 
in Fig. 5 below. 

The unprocessed date at box 46 are sent (as indicated by arrow 48) 
to a suitable display such as a CRT monitor as illustrated in Fig. 3 for graphic 
drsptay of tee waveform in real time (such as illustrated as grephic display element 
in Fig. 3). 

The data at box 46 are additionally analyzed using the user input 
parameters in accordance with the present invention. That is, the data may be 
used at box 52 to generate an ECG attractor (as set forth in greater detail in Fig. 9 
below) and then displayed as graphical display element 26 (Fig. 3). This graphical 
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T ■ reaMlme - ,rom * e orpu.se rate 

jal te visually represent me tempera, evotetton of oardiac dy„am,cs to mtfu. 

ZTTrT sys,ems exh,bit » a « * — 

d ^ fixed point or cross each other as the trajectories evolve over time, while 
pence* Rectories rollow a cydicai path. The date at box 46 may also be an! 
*zed te de.ee a Q RS event a, box 54 from which a RR time sequence or Z 
pulse sequence may be generated a, box 56 (as se, term in greater date,, in Fig. 6 

f M . oo, , ^ RR SeqUenC8 ( ° r inter -P ul! *> saquence) may be used for 
further HRV analysis, including determining a Chaotic Index at box 60 (as de- 
senbed ,n greater detail in Fig. 9 below,, and determining time domain parameters 
box 62 such as heart rate or pulse rate, RMSSD and SONN (as described in 
greater detail in Fig. 7) Ai, of these «me domain parameters may be dispiayed a. 
box 50 such as illustrated in Fig. 3. ^ 

TheRRtimesec l ue '" :e (orinter.pulsese<,uence)mayfurtherbeused 
a. box 64 to generate a heart rate (HR, or pulse rate (PR) time sequence or series 
(as described in greater detail in Fig. 7 betew). which may be sent (as indicated by 
arrow 66) for display as graphical display element 22 (see Fig. 3). Further, the HR 
time sequence or series may be used at box 70 for time-frequency (t-f) distribution 
analysis, including generating and displaying an intensity based color-mapped 
contour plot (as indicated by arrow 72) and generating (at box 74) and displaying 
(at box 50) SI. PI and Sl /P, indexes (as described in greater detail in Fig. 8 below) 

Fig. 5 illustrates the initiation of the analysis mode and data acquisi- 
tion. Specifically, the user first inputs the analysis moda at box 120. indicating 
whether operation is to use data being acquired at the time from a patient or 
operation is to use pre-acquired data. 
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If pre-acquired data is to be used, decision box 122 proceeds to list 
the availab.e files of pre-acquired data at box 1 24, and the user selects the desired 
file at box 126. If the user wishes to perform the analysis using a sampling fre- 
quency which is other than the default sampling frequency indicated at box 128 
(e.g., 500 Hz which detects the electrical signal of the heart 500 times per second) 
he may do so at box 130, in which case the sampling frequency may be changed 
to a different selected value at box 1 32. This may be required, for example, when 
the pre-acquired data of the se.ected file was acquired using a different sampling 
frequency than the default sampling frequency. Whatever sampling frequency is 
selected, the data from the selected file of pre-acquired data is then sequentially 
read at box 1 34 in time order. 

Alternatively, if analysis is to occur as the data is being acquired a 
determination may first be made by the user at box 140 as to whether or not the 
computer port receiving the data from the ECG data acquisition device is con- 
nected to the default port (e.g., a computer serial port such as previously de- 
scribed). In that case, if the user does not indicate at box 140 that a port different 
than the default port (e.g., COM Port 4) is to be used, then processing continues at 
box 142 with data acquisition occurring through the default port. If the user selects 
a different port, then the selected different port is set at box 144 to be recognized 
as receiving the data. Once the proper port for receiving data is set, the computer 
then begins to acquire data at box 146 from the ECG data acquisition device. 

As that data are acquired, whether from the computer file of pre- 
acquired data at a selected sampling frequency (at box 134) or from the ECG data 
acquisition device (at box 146), the data may at box 148 be displayed on the moni- 
tor 16 to showthe ECG waveform, which display may be updated periodically (e.g., 
every 0.1 seconds). 
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Of the da, /" SeqUenUa ' ^ aCqU ' red aCC0Klin9 '° *• ab0V ^ •»»— h» 
of the date then proceeds, including event detection 1 50 (F.g. 5, end nonlinear 

analyse 152 (Fig. 9, discussed further below). 

OR „ , . EV9m deteCtt ° n 38 musbatea "9- « 'nvolves determination of a 
QRS event ,n a sequential se« of date pointe in a time sedes of EC6 date, which 
may be characfenzed as Ma( t>. As is recognized by those stalled in me art, the 
ECG waveform of a standard heartbeat Is iilustrated in Fig. 1, with the standard 
peaks ,n that waveform having the conventional designations P, Q, R s and T 
Detection of a QRS event is the detecuon of an ECG waveform in the form o, 

ZZso 7T EnBelse ' W - AH - and Zee,enbe ^ c - : A s,n °* *« 

or QRS-DetecUon and Feature BrfraCon. Computers in Camlotogy 6, 37^.2 
1979 teaches QRS event detection, and the full disclosure thereof is hereby incor^ 
porated by reference herein. 

Initial filtering of the data first occurs. For example, a differentiator 
with a 62.5 Hz notch filter is applied at box 220, where: 

Y0(t) = data(t) - data(t-4) 
Such a differentiator filters out power line noise conventionally found at around 62 5 
Hz, as is explained in Friesen. G.M., Jannett, I.e., Jadal.ah, M.A., Yates S L 
Quint, S.R., and Nagle, H.T.: A Comparison of the Noise Sensitivity of Nine QRS 
Detection Algorithms. IEEE Transactions on Biomedical Engineering, BME-37 (1) 
PP85-98, 1990, the complete disclosure of which is hereby incorporated by refer- 
ence. In addition to filtering out power line noise, a low pass filter may also be 
applied at box 222 to filter out high frequency noise, where: 

Y1 (t) = Y0(t) + 4*Y0(t-1 ) + 6*Y0(t-2) + 4*Y0(t-3) +Y0(t-4) 
Such a suitable filter is also explained, for example, in Friesen et al 
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After the data have been suitably filtered, it is determined at box 224 
whe^ 

If t has not, then no R peak of a QRS event has occurred, in which case it is deter- 
ged at box 226 whether or not all of the current ECG data have been analyzed 
If there are more data for analysis, processing begins again at box 224 to deter- 
m«ne whether Y1 (t+1 ) has exceeded the threshold. ^ 

As the data Po^ts are sequentially processed according to the above 
,f the slope of the data points are found to exceed the threshold at box 224 it is 
then determined at box 230 whether or not crossing of the threshold occurred 
w,thm the search region. If it did not, then the high slope is interpreted at box 232 
as being a base.ine shift, such as can occur from, for example, breathing, move- 
ment of the patient, or a change in the contact of the electrodes with the patient's 
•fan. If it is determined at box 240 that the threshold crossings of the data points 
w,th,n the search region meet the established QRS event criteria (/.e., essentially 
■ndicates that the peak in the waveform is high enough to be an R peak) then the 
event is classified at box 244 as a QRS event. If it is determined at box 240 that 
the threshold crossings of data points within the search region do not meet the 
established QRS event criteria, then the event is classified as noise at box 246 

Once all of the current ECG data have been analyzed as determined 
at box 226 so that a QRS event has been detected, heart rate resampling and RR 
sequence generation (of the RR time series) proceeds at box 250. The RR time 
senes may then be used in non-linear analysis at box 152 as discussed further 
below in connection with Fig. 9. 

Heart rate resampling and RR sequence generation is illustrated in 
F,g. 7. At box 320, the location of R peaks as data points are identified (based on 
where QRS events were identified during event detection described above in 
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connection with Fig. 6). With the R peak data points identified, a sampling rate is 
chosen at box 322 for the heart rate signal (e.g., sampling frequency/100), such as 
is shown in Berger, R.D., Akselrod, S., Gordon, D., and Cohen, R.J., An Efficient 
Algorithm for Spectral Analysis of Heart Rate Variability. IEEE Transactions on 
5 Biomedical Engineering, BME-33 (9), pp900-904, 1986, the full disclosure of which 
is hereby incorporated by reference. The number of RR intervals (i.e., the interval 
from one R peak to the next R peak of a waveform) contained within a local window 
(e.g., 100 data points) of the heart rate signal is then calculated at box 326, and the 
instantaneous heart rate is calculated at box 328 as: 

10 HR = f * n/2, where f = sampling frequency, and n = number of RR intervals. 
The local window of data points is then slid forward (e.g., 100 points) at box 330 
and, if it is determined at box 332 that the end of the current ECG time series has 
not been reached, the number of RR intervals and the instantaneous heart rate is 
calculated for the new window of data points at boxes 328 and 330 as described. 

1 5 When it is determined at box 332 that the end of the time series has been reached, 
the heart rate time series may be displayed at box 336 and processing of a time- 
frequency distribution may proceed at box 340. 

In addition to the heart rate resampling in boxes 322-332, time do- 
main parameters may be simultaneously determined based on the location of the 

20 R peaks identified at box 320. The calculation of time domain parameters is usu- 
ally based on the NN interval, as opposed to the RR interval. Therefore, potential 
ectopic beats should be removed from the calculation. In this calculation therefore, 
ectopic beats may be identified and removed at box 360, for example by determin- 
ing that if the RR interval duration is more than 1 5% different than the previous RR 

25 interval, such a beat would be taken to be ectopic and then removed. That is, if 
abs[RR(i)/RR(i-1 ) - 1] > 0.1 5, then the beat at time V is taken as being ectopic and 
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.s removed. Checking for ectopic beats, and removing those identified ectopic 
beats at box 350 is continued until it is determined at box 352 that data with ectopic 
beats removed have been collected for a five minute period, at which point SDNN 
and RMSSD may be calculated at box 354 as follows: 

RMSSD = sqrt(sum(sqr(rr(n)-rr(n-1)) t N)/N) 
SDNN = stdev(n. N) 
After SDNN and RMSSD have been calculated, they are displayed at box 356. 

Determination of time-frequency distribution is shown in Fig 8 Joint 
time-frequency distributions may be used to depict the time-varying behavior of 
s.gnals of which the frequency content is of interest. Use of one of the Wigner-Ville 
family of time-frequency distributions make it possible to achieve an appropriate 
compromise between frequency resolution and time resolution. As illustrated in 
Fig. 8, a uniformly sampled HR time series is obtained at box 420 such as previ- 
ously described in connection with boxes 320-340 (Fig. 7). The HR time series is 
1 5 then converted at box 422 to an RR time series, where: ' 

RR(t) = 60,000/HR(t) 
The time-frequency distribution for the RR time series is then calculated at box 426 
using the kernel function which is empirically determined to be optimal, such as 
described in Pola, S., Macerata. A.. Emdin, M.. and Marchesi, C. Estmation of the 
Power Spectral Density in Nonstationary Cardiovascular Time Series: Assessing 
the Role of the Time-Frequency Representations (TFR). IEEE Transactions on 
Biomedical Engineering, Vol. 43, No. 1 , pp 46-49, the complete disclosure of which 

is. hereby incorporated by reference 

Jointt-f distribution analysis mathematically decomposes the RR time 
series into time-varying components of the frequency spectra. These frequency- 
domain calculations result in three main HRV spectral components: very low fre- 
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quency (VLF). low frequency (LF), and high frequency (HF). The LF component 
(0.04 to 0. 1 5 Hz) has been associated mainly with sympathetic activity while the HF 
component (0.15 to 0.40 Hz) has been correlated with parasympathetic activity. 
There is a constant interplay between these autonomic stimuli to influence HR. The 
5 resulting sympatho-vagal balance can be quantified by using the ratio of LF to HF 
spectral power. In this context, analysis using frequency methods has been found 
to be a better predictor of physiologic changes than time-domain methods. 

An SI index, PI index and SI/PI ratio may be calculated (and dis- 
played) at boxes 430, 432 and 436. The SI index is the spectral power in the 0.04 

10 Hz to 0. 1 5 Hz low frequency range of the t-f distribution integrated over the entire 
time duration of the t-f distribution displayed on the computer screen, the PI index 
is the spectral power in the 0.15 Hz to 0.40 Hz high frequency range of the t-f 
distribution integrated over the entire time duration of the t-f distribution displayed 
on the computer screen, and the SI/PI ratio is the ratio of SI spectral power to PI 

15 spectral power. The SI/PI ratio is a quantification of the above mentioned 
sympathovagal balance. The SI index, PI index and SI/PI ratio correspond to the 
moment by moment predominantly sympathetic tone, parasympathetic tone, and 
sympatho-vagal balance, respectively. 

The time-frequency distribution for the RR time series calculated at 

20 box 426 may also be color mapped according to spectral power intensities at box 
450, and the intensity-mapped color display representing that distribution may be 
displayed at box 452. More specifically, the color mapping consists of converting 
the time-frequency (t-f) distribution values to color-coded intensity maps which have 
been found to visually illustrate certain data conditions which a physician may find 

25 useful. This color mapping may be accomplished by determining the maximum 
value (global max) of the t-f distribution for the entire time and frequency range for 
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the RR interval sequence being analyzed. The frequency range may be fixed to 
limits corresponding to the ranges used for computing SI and PI, namely: 

Hf_max = 0.4 Hz 

Hf_min = 0.15Hz 

5 Lf_max = 0.15 Hz 

Lfjnin = 0.04 Hz 
hfreq = 0.5 Hz 

For each time slice to be displayed, the local minimum and maximum values of the 
t-f distribution are evaluated over the frequency range from Lfjnin to hfreq (e.g., 

10 from 0.04 Hz to 0.5 Hz), with the local minimum value being defined as the single 
lowest power value in that frequency range and the local maximum value being the 
data point with the highest power value that is greater than the power of its two 
preceding and two succeeding neighbors in the range. A weighted average time- 
slice maximum (M WA ) is then calculated: 

1 5 M WA = 0.9 local max + 0.1 Global max 

For each time slice, the color value for each point in the t-f distribution is calculated 
as follows: 

Color value = (Power value - local min) / (M WA - local min) 
The color value is then mapped using a table of, for example, 256 rainbow colors 
20 ranging from black to white. The following is an example of a table which may be 
suitably used to map the color value in RGB (Red, Green, Blue) format: 
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Graphic display element 24 illustrates an intensity-mapped time-frequency distribu- 
tion color contour plot, with its time axis (the horizontal axis) shared with the time 
axis of display element 22. 

Non-linear analysis may also be performed using the RR time series 
(from box 250, Fig. 6) and data waveform (from box 148, Fig. 5) as illustrated in 
Fig. 9. 

Specifically, one analysis which may be performed is to use the data 
point time series for the waveform, whether received from the ECG device or a pre- 
acquired data file (from box 148. Fig. 5), to generate at box 522 a XY scatter plot 
conventionally known as an "ECG attractor". As is known to those skilled in the art, 
the original ECG time series (i.e., data(t)) is used as the Y-coordinate and its time^ 
embedded equivalent time series (i.e., data(t-tau)) is used as the X-coordinate, 
where tau equals two ECG sample intervals. That is, a delay filter is used to gener- 
ate the nth dimension of data from the (n-1 ) dimension. In the illustrated example, 
a static delay of two sample intervals is used to generate the second dimension, 
although it should be understood thatthis can be extended to more dimensions and 
different delays. The XY scatter plot of the ECG attractor is displayed at box 524. 

In another similar non-linear analysis (using the RR time series from 
box 250, Fig. 6, whereas the graphical representation of the ECG attractor of box 
522 uses the raw ECG waveform from box 148, Fig. 5), a two-dimensional XYtime 
series characterized as an RR attractor is generated at box 530 using the RR time 
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series as the Y-cbordinate and its time embedded equivalent time series as the X- 
coordinate, with a time delay equal to two RR sample intervals. The nearest neigh- 
bor and its separation from the initial point is then determined at box 532. If the 
spatial separation is determined at box 534 to be greater than a selected threshold, 
5 then Gram-Schmidt reorthonormalization is performed at box 536 on the vector 
defined by the two points and then the step of box 532 is repeated until the separa- 
tion does not exceed the threshold for renormalization (as determined at box 536), 
at which point the principal axis vector can be obtained at box 538. The principal 
axis vector may then be used at box 540 to estimate the largest Lyapunov expo- 
1 0 nent, conventionally known as the Chaotic Index, which is a measure of the degree 
of chaos. Algorithms for calculating Chaotic Index are known in the art, such as 
. shown in Wolf MM, Varigos GA, Hunt D, Sloman, JG. Sinus arrhythmia in acute 
myocardial infarction. Med. J. Aust, 2 pp 52-53, 1978, the complete disclosure of 
which is hereby incorporated by reference. 
15 The Chaotic Index may then be displayed at box 542. Determining 

and periodically monitoring the Lyapunov exponent for a physiological system over 
an extended length of time could reveal additional trends towards less or more 
chaotic behavior which may be indicative of a progressive disease requiring phar- 
maceutical or therapeutic intervention or an adjustment of a treatment regimen. 
20 In accordance with the present invention, the above data analysis 

may be performed in a dynamic manner by refreshing the analysis in real-time 
(where real-time is used herein as referring not only to analysis occurring while 
external data are being received but also to dynamic analysis of pre-acquired data 
as those data are played back over a time period essentially corresponding to the 
25 time period for which the data were previously acquired). 
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More specifically, as indicated in Fig. 4, the ECG or pulse waveform 
may be continuously updated to display a waveform at box 50 based on, as an 
example, a batch of 4,096 data points received at box 46. At a 500 Hz sampling 
frequency (comprising 500 samples per second), 4,096 data points comprise 
5 approximately eight seconds of data. Those same 4,096 data points are also XY 
plotted at box 52 and displayed at box 50 as the ECG attractor. 

Those same 4,096 data points are also processed for event detection 
(see Fig. 6) at box 54, with a continuous output of detected QRS events used to 
generate a continuous RR sequence (RR time series) at box 56 (see box 250, Fig. 

1 0 6). With event detection being accomplished within the time frame of the batch of 
data points being processed (e.g., approximately eight seconds for 4,096 data 
points at a 500 Hz sampling frequency), receipt and event detection for the next 
batch of 4,096 data points may be accomplished in real-time (i.e., in keeping with 
the time element of the data points). 

1 5 While the event detection may be accomplished with a batch of data 

points, each detected RR interval may nevertheless be output to box 56 (generat- 
ing the RR time series) as it is detected {i.e., before event detection is completed 
for all 4,096 data points). Therefore, it should be appreciated that at this point in 
the continuous processing after box 56 in Fig. 4, analysis will occur using RR 

20 intervals as "points" rather than raw data points as used in box 54. 

Specifically, the Chaotic Index may be determined at box 60 for every 
128 RR intervals as determined at box 56. Therefore, once the first 128 RR inter- 
vals have been determined, the Chaotic Index will be calculated (see box 540, Fig. 
9) and displayed until another group of 128 RR intervals have been determined, at 

25 which point a new Chaotic Index will be similarly calculated and the displayed 
Chaotic Index will change to the newly calculated Chaotic Index. 



10 



15 



20 



25 



Provisional Application 
01819P0080US 



-24- 



a. bo* 62 are also calculated using RR intervals as date points These 1^ 
-^^atedanddMayedwn.n^n,^^^^^ 
been accumulated, and then may be recalculated and displayed KeS^Z 
ume a new RR inten,a, date point is received. Thus, the SDNN param^e 

■me seamen, o, elecfcocandiogram (or pulse pleteysmograph) date after puTl 
abnormal RR intervals (or inter-pulse intervals) ana removed (a. box 3o0 n Ro 7 T 
and the RMSSD parameter i. ih« , . 9 ' 7) ' 

successive RR , , of the difference between 

success,ve RR intervals (or ,nter-pulse interva.s) fram .he same 5 minute segment 
of electrocardiogram (or pulse plethyamograph) data. 

which te ProCeSS ' ng bey °" d b ° X 64 USSS HR se<)uenoe < HR «™ -nee) 
wh,ch s a conversion of the RR time series to a time series having a Uniterm 

interval see box420. Rg. 8, Thosedate may Ihen oe displayed (peranToJte 
display 256 uniform interval points accumulated in forty poin, belches. The date 
may s,m„ariy be used in such 256 point groupings forlime-ftequenoy (,-f, dislribu! 
on analyst* from which tee SI, P, and Sl/P, indexes may be calcu lated a, 
(see boxes 430^36 of Rg. 8) and men displayed and tee ,-, distention spectra, 
power .ntensiliea color mapped and displayed (boxes 45(M52 of Rg 8 ) Color 
mapprng and index calculate*, may be performed using a moving window o, 256 
HR ,nte,val po,nts (from box 64). updated every 64 points. That is. when 256 HR 
rnterval pornts are firs, raceived. mapping and index calculation will occur with tee 
rasute delayed until an additional 64 HR Inte™, points are received, a. whfch 
pom, tee flrs, 64 HR interval points in the prior batch will be dmpped and .hen 
mapprng and index ca,cu.a.ion will again be done using tee las. 192 HR leZ 
perns from .he prior batch and the new 64 HR interval points 
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It should be recognized that the present invention is not limited to the 
above details relating to suitable processing of data points, including the particular 
numbers of points used in individual calculations. However, it should be appreci- 
ated that the above described manner of processing the received heart beat data 
5 has been found to be suitable in providing the desired real-time analysis, with the 
attendant advantages to physician knowledge and patient care. It should also be 
appreciated that instead of patient treatment, the method and system of the pres- 
ent invention could be used for alternative purposes, such as relaxation training. 

It should also be appreciated that much of the above may be accom- 

10 plished using suitable software performing the described processing and display. 
Visual C++ is one suitable programming language which may be used, as illus- 
trated by Figs. 10A-10C which sets forth a Visual C++ programming language 
algorithm for detecting the peaks in the R wave and computing the R-R interval (for 
QRS event detection, as shown in Fig. 6). 

15 F '9s. 11-12 illustrate screenshots displaying the data and analysis 

processed according to the above description for abnormal conditions, with Fig. 1 1 
illustrating an example HRV analysis on pre-acquired representative data for a 
Chronic Heart Failure (CHF) electrocardiogram and Fig. 12 illustrating an example 
HRV analysis on pre-acquired data representative of an epileptic seizure episode 

20 electrocardiogram. As can be seen in Fig. 1 1 , the HRV analysis for the example 
Chronic Heart Failure (CHF) electrocardiogram results in a display of a heart rate 
of 65 beats per minute, an SDNN of 36 ms, an RMSSD of 27 ms, an SI/PI ratio of 
0.4, and a Chaotic Index of 0.04. In Fig.. 12, by contrast, where the input data 
involve an epileptic seizure episode, the heart rate is 60 beats per minute, the 

25 SDNN is 238 ms, the RMSSD is 30 ms, the SI/PI ratio is 1 .4, and the Chaotic Index 
is 0.60. The differences in the displayed information between the normal condition 
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of Fig. 3 and the different abnormal conditions of Figs. 11-12 (the various displayed 
indices, as well as the displayed plots [e.g., electrocardiogram, heart rate, intensity- 
mapped time-frequency distribution color contour plot, electrocardiogram attractor]) 
provide an important new tool for the development of a detailed understanding of 
the dynamic mechanisms underlying the conditions represented, and may provide 
distinct and valuable data to a treating physician who, when provided in real-time 
as the patient undergoes the abnormal condition, can be assured of having the 
most up to date information for evaluation as he evaluates possible treatments. 

The system and method of the present invention provide a new tool 
for real-time automated analysis of heart rate variability and its adjuncts, that 
combine the non-stationary analysis capability for evaluating cardiac signal histo- 
ries with the predictive capability of non-linear analysis to better monitor and cate- 
gorize autonomic regulation of cardiac function. 

Still other aspects, objects, and advantages of the present invention 
can be obtained from a study of the specification and the drawings. It should be 
understood, however, that the present invention could be used in alternate forms 
where less than all of the objects and advantages of the present invention and 
preferred embodiment as described above would be obtained. 
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CLAIMS 

1. A method of monitoring variability of heart activity occurring 
2 during a time period, comprising the steps of: 

sequentially receiving data points of heart activity data over a period of time 
4 corresponding to the time period of the heart activity; 

evaluating said data points as sequentially received to determine QRS 
6 events; 

outputting said QRS events to a processor as they are sequentially deter- 
8 mined; 

processing said output QRS events as they are output to periodically deter- 
1 0 mine heart rate variability information, wherein said heart rate vari- 

ability information is based on a selected number of output QRS 
12 events; 

periodically redetermining said heart rate variability information using at 
14 least some subsequently output QRS events; 

during said period of time corresponding to the time period of the heart 
1 6 activity, displaying the most recently determined heart rate variability 

information. 
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2. A method of monitoring variability of heart activity occurring 
2 during a time period, comprising the steps of: 

sequentially receiving data points of heart activity data over a period of time 
4 corresponding to the time period of the heart activity; 

evaluating said data points as sequentially received to determine QRS 
6 events; 

outputting said QRS events to a processor as they are sequentially deter- 
8 mined; 

processing said output QRS events as they are output to repeatedly defer- 
. 1 0 mine a Chaotic Index of a selected group of determined QRS events, 

wherein said selected group of determined QRS events comprises X 
1 2 QRS events with said Chaotic Index being determined for sequential 

groups of X QRS events; 
14 during said period of time corresponding to the time period of the heart 

activity, displaying the most recently determined Chaotic Index, 



3. The method of claim 2, wherein X is 128. 
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4. A method of monitoring variability of heart activity occurring 
2 during a time period, comprising the steps of: 

sequentially receiving data points of heart activity data over a period of time 
4 corresponding to the time period of the heart activity; 

evaluating said data points as sequentially received to determine QRS 
6 events; 



outputting said QRS events to a processor as they are sequentially deter- 



processing said output QRS events as they are output to repeatedly deter- 
mine time domain parameters, wherein said time domain parameters 
are first determined after a selected portion of said time period and 
repeatedly updated thereafter additionally using subsequent second 
selected groups of QRS events; 

during said period of time corresponding to the time period of the heart 
activity, displaying the most recently determined time domain param- 



8 



mined; 



16 



eters. 



5. 



The method of claim 4, wherein the determined time domain 



2 



parameters are SDNN and RMSSD. 
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2 ' h,,h 6 " A meth ° d ° f m ° nitorin 9 variabj "ty of heart activity occurring 

2 dunng a time period, comprising the steps of: 

sequentially receiving data points of heart activity data over a period of time 

corresponding to the time period of the heart activity 
evaluating said data points as sequentially received to determine QRS 
° events; 

outputs said QRS events to a processor as lhey are ^ 

minor!* 



8 mined; 



processing a selected number of output QRS events to determine a heart 
rate time series having a uniform interval for said selected number of 
most recently output QRS events, wherein said determined heart rate 
time series is updated after an additional Y number of QRS events 
are output using the most recently output selected number of QRS 



14 events; 



processing said determined heart rate time series to determine a time-fre- 
quency distribution for the most recently updated heart rate time 
series; 

displaying the most recently determined time-frequency distribution. 

7. The method of claim 6, wherein said selected number is 256. 

8. The method of claim 7, wherein Y is 40. 

9. The method of claim 6, wherein said displaying step display 
uses intensity-mapped colors. 
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10. A method of monitoring variability of heart activity occurring 
2 during a time period, comprising the steps of: 

sequentially receiving data points of heart activity data over a period of time 
4 corresponding to the time period of the heart activity; 

evaluating said data points as sequentially received to determine QRS 
6 events; 

outputting said QRS events to a processor as they are sequentially deter- 
8 mined; 

processing a selected number of output QRS events to determine a heart 
1 0 rate time series having a uniform interval for said selected number of 

most recently output QRS events, wherein said determined heart rate 
12 time series is updated after an additional Y number of QRS events 

are output using the most recently output selected number of QRS 
14 events; 

processing said determined heart rate time series to determine a time-fre- 
16 quency distribution for the most recently updated heart rate time 

series; 

18 processing the most recently determined time-frequency distribution to 

determine its spectral power in a low frequency range and its spectral 
20 power in a high frequency range of the t-f distribution; 

displaying the most recently determined spectral power in the low frequency 
22 range and the spectral power in the high frequency range. 

1 1 . The method of claim 10, wherein said selected number is 256. 



12. The method of claim 1 1 , wherein Y is 40. 
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13. The method of claim 10, displaying the ratio of the most re- 
cently determined spectral power in the low frequency range to the most recently 
determined spectral power in the high frequency range. 

14. The method of claim 10. wherein the low frequency range is 
0.04 Hz to 0.15 Hz. 

1 5. The method of claim 1 0, wherein the high frequency range is 
0.15Hzto0.4Hz. 

16. A method of monitoring variability of heart activity occurring 
during a time period, comprising the steps of: 

sequentially receiving data points of heart activity data over a period of time 

corresponding to the time period of the heart activity; 
evaluating said data points as sequentially received to determine QRS 
8 events; 

outputting said QRS events to a processor as they are sequentially deter- 
10 mined; 

processing said output QRS events as they are output to repeatedly deter- 
mine a Chaotic Index of a selected group of determined QRS events, 
wherein said selected group of determined QRS events comprises X 
QRS events with said Chaotic Index being determined for sequential 
groups of X QRS. events; 
during said period of time corresponding to the time period of the heart 
activity, displaying the most recently determined Chaotic Index; 
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processing said output QRS events as they are output to repeatedly deter- 
mine time domain parameters, wherein said time domain parameters 
are first determined after a selected portion of said time period and 
repeatedly updated thereafter additionally using subsequent second 
selected groups of QRS events; 

during said period of time corresponding to the time period of the heart 
activity, displaying the most recently determined time domain param- 
eters; 

processing a selected number of output QRS events to determine a heart 
rate time series having a uniform interval for said selected number of 
most recently output QRS events, wherein said determined heart rate 
time series is updated after an additional Y number of QRS events 
are output using the most recently output selected number of QRS 
events; 

processing said determined heart rate time series to determine a time-fre- 
quency distribution for the most recently updated heart rate time 
series; 

displaying the most recently determined time-frequency distribution; 

processing .the most recently determined time-frequency distribution to 
determine its spectral power in a low frequency range and its spectral 
power in a high frequency range of the t-f distribution; 

displaying the most recently determined spectral power in the low frequency 
range and the spectral power in the high frequency range. 

1 7. The method of claim 1 6, wherein said selected number is 256. 



2. 



8 
10 



16 
18 



Provisional Application 
01819P0080US 



-34- 

1 8. The method of claim 1 7, wherein Y is 40. 

19. The method of claim 16, wherein said data points of heart 
activity data are received during the heart activity. 

20. The method of claim 16. wherein said data points of heart 
activrty data are received from a pre-acquired file of data points of the heart activity. 

21 . A system for monitoring variability of heart activity occurring 
2 during a time period, comprising: 

a heart activity data acquisition device adapted to acquire sequential data 
4 points of heart activity of a patient; 

memory adapted to store sequential data points of heart activity in pre-ac- 
6 quired data files; 

a user input for selecting between said acquisition device and a selected 

pre-acquired data file as a data source; 
a processor adapted to 

from said selected data source, sequentially receive data points of 
heart activity data over a period of time corresponding to the 
12 time period of the heart activity, 

determine QRS events from said data points as sequentially re- 
14 ceived, 



output said QRS events as they are sequentially determined, 
repeatedly determine a Chaotic Index of a selected group of deter- 
mined QRS events as they are output, wherein said selected 
group of determined QRS events comprises X QRS events 
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with said Chaotic Index being determined for sequential 
20 groups of X QRS events, 

repeatedly determine time domain parameters from said QRS events 
22 as they are output, wherein said time domain parameters are 

first determined after, a selected portion of said time period 
24 and repeatedly updated thereafter additionally using subse- 

quent second selected groups of QRS events, 
26 determine a heart rate time series having a uniform interval for a 

selected number of most recently output QRS events, wherein 
28 said determined heart rate time series is updated after an 

additional Y number of QRS events are output using the most 
30 recently output selected number of QRS events, 

determine a time-frequency distribution forthe most recently updated 
32 heart rate time series, and 

for the most recently determined time-frequency distribution, deter- 
34 mine spectral power in a low frequency range and its spectral 

power in a high frequency range; 
36 a display continuously updated during said period of time corresponding to 

the time period of the heart activity to display the most recently deter- 
38 mined Chaotic Index, the most recently determined time domain 

parameters, the most recently determined time-frequency distribution, 
40 the most recently determined spectral power in the low frequency 

range, and the most recently determined spectral power in the high. 
42 frequency range. 



22. The system of claim 21 , wherein X is 128. 
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23. The system of claim 21 , wherein said selected number is 256. 

24. The system of claim 23, wherein Y is 40. 



25. The system of claim 21 , wherein processor includes function 
2 calls for Fast Fourier Transforms and Inverse Fast Fourier Transforms. 
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HRV/PRV Real-time Analysis System 
Software and Data Flowchart 



User Input 
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ExtamaJ signal 



44 



-A. 



Pto-tq cquItc d data 



Unprocessed EKG or Pulse 
Pleftysmograph Input data 



4096 pts 



64- 



Event DatecSon 



RR or Inter-pulse sequence 



KR or PR sequence 



256 pts. accumulated in 
40 pt batches 



lb 



W Distribution 



4096 
pts 



EKG or Pulse Rethysmogmph time 



•attractor 



EKG or Pulse plathysmoarapti waveform 



4096 pts 



<3 



60 
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CnaoSc index 



Updated once every 128 RR pts 



Time domain parameters- 
HRorPR. RMSSD.SDNN 



Walts for first 5 min of RR data for initial 
calculation then updates every point 

66 



HRorPR Time Series 



256 pts accumulated in 40 
- pt batches 

intensity based color-mapped contain- plot 



256 pts updated once every 64 t r f pts 



SI, Pt end SI/PI Indexes 



Updated once every 64 t-f pts 
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USER INPUT AND ECG ACQUISITION 



Pre-acquired 
data files 



Yes 




IZO 




User Input* 
t Filename selection 




/ 



Use default 500Hz 
sampling 
frequency 



Yes 



Change sampling 
frequency to user 
selected value 




iM- 



Change COM Port 
to user selected 
value 



5 default COM 
Port 4 



Read data front 
selected pre- 
acquired file 



7 



Start acquiring 
data from external 
EKG device 



/Diplay i 
/ EKG\ 



and refresh 
waveform 
V after every 0.1- 





Event 
Detection 



QRS EVENT DETECTION 




EKG Time series: data(t) 



Apply differentiator with 
62.5H2 notch filter (see 
Friesen et al. reference) 

Y0(t) = data(t)-data(t-4) 



•no- 



Apply low pass filter (see 
reference ?) 

Y1(t) = Y0(t) + 4*Y0(t-1) 
+6*Y0(t-2) + 4*Y0(t-3) + 
- YO(M) 




HEART RATE RESAMPLING ANriF°°^ s ^ '^ ^3 -'Q3KB« 
SEQUENCE GENERATION ALGORITHM 




2So 
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Identify location of R 
peaks as data points 
where QRS events 
were identified during 
event detection 



Identify and remove ectopic 

beats using the criteria 
8bs[RR(iyRR(i-1) - ij > 0.15 



Choose sampling rate 
for heart rate signal as 
sampling frequency/ 
100 

(see Berger et a!. 
reference) 



3Zfe 



35^ 



Calculate number of 

RR intervals 
contained within a 
100 data point local 
window of the heart 
rate signal . 




Calculate 
instantaneous heart 

rate as HR 
0.5'samplingfreq. 
•Number of RR 
intervals 



Slide the 100 data 
point local window 
100 heart rate data 
points forward 



Yes 



Calculate RMSSD 
and SDNN 



— 33o 



340 




-354- 





TIME-FREQUENCY DISTRIBUTION 




HR time series from Rate 
Resampling algorithm output 



Convert HR time series to 
RR time series: 



RR(t) = 600Q07HR(t). 



Calculate time-frequency 
distribution for RR time 
series using exponential 
kernel function (see 



Integrate the t-f distribution 
spectral power m the 0.04Hz 
to 0.15Hz range and display 
value as SI index • 



Integrate the t-f distribution 
spectral power in the 0.1 5Hz 
to 0.40Hz range and display 
value as PI Index 



Compute and display the 
ratio of the SI and PI indexes 



Color map the time-frequency 
distribution spectral power . 



46c 



Display the intensity- 
mapped color display 
representing the time- 
frequency distribution 
< for the RR time series 
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NON-LINEAR ANALYSIS 



EKG device or 
p re-acquired data 
waveform 



data(t) 




Generate a XY scatter plot 
("ECG attractor") by using 
original EKG time series i.e. 
data(t) as the Y-coordinate 

and Its time-embedded 
equivalent time series data(t- 

tau) as the X- coordinate. 
The time delay tau = 2 EKG 
sample intervals 



•S9- ' iiA 



. . Display 
^ "ECG Attractor" 



RR time series 



531. 



I 



Generate a 2-dimensional 

XY time series ("RR 
attractor*!) by using original 
RR time series as the Y- 
coordinate and its time- 
embedded equivalent time 
series as the X- coordinate. 
The time delay = 2 RR 
sample Intervals 



For each point in the "RR 
attractor" time series 
. determine its closest 
neighbor and their 
separation 
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Perform Gram- 
Schmidt 
reorthonormaHzation 



55o 
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Obtain principal axis 
vector 
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Estimate Largest 
Lyapunov Exponent 
"Chaotic Index" • 
(see Wolf et al. 
reference ?) 



/ Oh 



Display 
\ "Chaotic Index" 



Figure 10A 

EventStream::DetectEvents() // dfl version 

mt*yO; 
int*y1; 

bool eventoocunred=1alse; 

yO=(int *)marioc(MAXDATASI2E*sizeof(int));. 
y1=(int nmalloctMAXDATASIZ^sizeoflint)); 



fbr(int N4; i<MAXDATASIZE;i++) { 
y0[i]=datap]-datap-4]; 



for (I=8;KMAXDATASIZE;i++) { 

yin=y0[Q+4*y0pw1]+6*y0n-2]+4*yOP-3J+yOti-4]; * - , 

y1I0l=data[0>prevdata[4]+4*(prevdatan^revdata[3])+ 

6*(prevdata[6]-prevdataI2])+4*(prevdata[5]-prevdata[1])+ 
prevdata[4Kprevdatap)]; 

y1 [1]=data[1]-prevdata[5]+4*(data[0]-prevdata[4])+ 
. 6*(prevdataI^prevdata[3])+4*(prBVdata[6h)revdata[2])+ 
prevdata[5hprevdata[1 ]; 

y1[2f=datap]-prevdata[6l+4*(data[1]-prevdatat5])+ 

6*(data[01-prevdataI4])+4*(prBvdata[7]-prevdata[3])+ 
prevdata[6]-prevdata[2J; 

y1 [3]=data[3J-prevdata[7]+4*(data[2]-prevdata[6J>+ 
6*(data[1]-prevdataI5])+4*(data[0]-prevdatat4])+ 
prevdata[7]-prevdata[3]; 

y1[4J=dala[4]-data[0]+4*(data[3J-prevdata[7])+ 

6*(data^prevdata[6])+4*{data[1]-prevdata[5])+ 
data[0]-prevdata[4]; 

y1[5J=data[5Hata[1]+4*(data[4]-data[0])+ 

6*(data[3J«prevdata[71H4*(dataP]H3revdata[6])+ 
data[1J-prevdata[5]; 

y1 [6]=data[6hdata[2]+4*(data[5]-dataL1])+ 

6*(data[4J-data[0])+4*(data[3].prevdatar7])+ 
data[2>prevdata[6]; 

y 1 r71=datap>data[31+4*(data[q-data[2])+ 

6*(data[6>data[1])+4*(data[4l-prevdata[0])+ 
data[3J-prevdata[7]; 

if (detectcalls=0){ 

iastevent=0; ■ " " 

prevlastevent=0; 

curevent=0; 

} 



Figure tOB 

eventcounter=0; 

for (Int j=0J<MAXDATASI2E;j++) 
eventanray[j]=0; 

fbr(r=0;i<58;i++){ 
Sf(prevy1p]>200){ 

forQ=0;j<58;j++){ 
if(i+j>57){ 

if(yip+f58]<-200) 



} 

else 



eventoccurred=true; 



eventoccurred=false; 
if (eventoccurred) 
break; 



lf(prevyin+jl<-200) 

eventoccurred=true; 



}//end fbrj 

} 



eventoccurred=faise; 
if (eventoccurred) . . 

break; 



eventoccurred=false; 
if (eventoccurred) { 

if ((detectcaI!s*MAXDATASIZE+i-58-curevent)>100) { 

prevlastevent=lastevent; 

lasteventFCurevent; 

curevent=detectcalls*MAXDATASI2E*H58; 
. if (curevenWastevent>10Q0) { 
lastevent=curevent; 
previastevent=curevent; 

} - 
>{ 

rrarray[rrcounter3=curevent-iastevent; 
rrcounter++; 



} 

}// end fori 



for (i=0;i<MAXDATASIZE-58;i++) { 

if(y1[i]>200){ 

for(j=0;j<58;]++){ 

if(y1D+fl<-200){ 

eventoccurred=true; 

} 



eventoccurred=false; 
if (eventoccurred) 
break; 

}//endforj 

else 

eventoccurred=fafse; 



Figure 10C 



if(eventoccurred){ 

ff((detectcaIIs*MAXDATASIZE+i-curevent)>100) {// High HR reject 

prevlastevent=lastevent; 
lastevent=curevent; 

curevent=detectcat!s*MAXDATAS(ZE+i; 

if (curevent-lastevent>1 000) { // low HR reject * 
Iastevent=curevent; 
prevlastevent=lastevent, 

} 

efse { 

rrarraylrrcounterl=curevenWastevent; 
rrcounter++; 

;} . ' . 

while (deteCt^lls*MAXDATAS!2E+eventcounter<curevent-50 
&& prevlastevent<lastevent-1Q0 

. && lasteve nKcurevent-1 00) { // exception for first run and skips 



prevfastevent); 



lastevent); 



if (detectcails*MAXDATASIZE+eventcounter<iastevent-50){ 

eventarray[eventcounter]=60.0Tsamplerate/(lastevent- 

) 

.else if (detectcaJls*MAXDATASfZE+eventcounter<lastevent+50K 
int portiortinlast=lastevent-(detectcaHs*MAXDATASI2E+ 

eventcounter); 
float numberofrrinterva!s= 

(50.0f*(float)portioninlast)/ 
((float)lastevent-(float)prev!astevent)+ 
(50.0f-{float)portioninlast)/ 
((float)curevent-(float)lastevent); * 
eventarray[eventcounter]=numberofrrfntervals/ 
(1 00.0f/(float)samplerater60.0f, 

.} 

else ff (detectc«IIs*MAXDATASIZE+eventc6unter<curevent-50) { 
eventarrayteventcounter}=60.0rsamplerate/(curevent-- 

> 

. eventcounter+=100; 
}// end while 



}//endff 
}//endif 
} // end for 

;for(i=0;i<8;l^+) 

prevdatap]=data[MAXDATA§fZE-8+l]; 

for(i=0;i<58;i++) 
. prevyip]=y1IMAXDATAS>ZE-58+iJ: 

detectcalls++;« 
free(y0); 
free(y1); 
retum(O); 
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